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Preface 


The  purpose  of  this  research  is  to  investigate  two 
methods  of  spectral  estimation  for  typical  radar  type 
signals.  The  two  methods  investigated  are  Blackman-Tukey 
and  Burg.  A  comparison  of  the  two  methods  is  made  based  on 
their  statistical  properties  and  computational  efficiency. 

Valuable  guidance  was  provided  by  Maj  Glenn  Prescott  of 
the  Air  Force  Institute  of  Technology  (AFIT)  during  the 
course  of  this  work.  His  help  is  gratefully  appreciated. 

My  deepest  gratitude  is  expressed  to  my  dear  wife,  Alfreds, 
for  her  constant  support  and  encouragement.  Her  help  in 
preparing  the  manuscript  has  enable  me  to  accomplish  a  major 
goal  in  my  life.  Another  major  goal  that  we  both 
accomplished  during  my  stay  at  AFIT  was  the  birth  of  our 
beautiful  daughter  Ashlee.  A  special  thanks  is  given  to  our 
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Abstract 


The  purpose  of  this  study  was  to  examine  the 
Blackman-Tukey  (BT)  and  Burg  methods  of  spectral  estimation 
for  typical  electronic  warfare  received  signals.  Such 
signals  are  generally  short  in  duration,  resulting  in  short 
data  records.  The  BT  method  is  a  conventional  spectral 
estimation  scheme  and  is  based  on  computing  the  discrete 
Fourier  transform  of  the  autocorrelation  sequence  (ASC) 
derived  from  the  data  record.  An  inherent  problem  of  this 
approach  is  that  of  data  windowing.  Data  windowing  may 
result  in  poor  frequency  resolution,  particularly  for  short 
data  records . 

The  Burg  method  of  spectral  estimation,  a  modern 
approach,  is  capable  of  providing  relatively  good  frequency 
resolution  for  short  data  records.  However,  this  method 
requires  sufficient  input  signal-to-no ise  ratio  (SNR).  The 
idea  here  is  to  extend  the  ACS  by  extrapolation  (or 
prediction)  rather  than  windowing  the  data. 

The  Burg  method  was  found  to  yield  far  superior 
performance  for  data  records  consisting  of  64  data  samples. 
Note,  however,  that  a  minimum  SNR  of  15  dB  was  assumed. 
Using  this  method  a  "smart"  routine  was  developed  that 
automatically  determines  the  actual  frequency  components  of 
the  data  record . 


SPECTRAL  ANALYSIS 


OF  SHORT  DATA  RECORDS 

I  .  Overview 


Introduction 

The  motivation  for  this  study  is  the  need  for  improved 
detection  of  hostile  signals  (i.e.,  signals  that  are 
radiated  by  hostile  emitters).  An  emitter  is  any  signal 
source  (i.e.,  a  radar,  a  jammer,  etc).  The  signal  radiated 
by  the  emitter  is  made  up  of  many  measurable 
characteristics.  Wiley  (33:8-12)  provides  a  summary  of  the 
measurable  signal  characteristics  and  the  emitter 
capabilities  inferred  from  these  characteristics.  The  most 
common  measurable  signal  characteristics  are: 

1.  Pulse  repetition  frequency  ( PRF ) 

2.  Transmitted  signal  power 

3.  Transmitted  beamwidth 

4.  Operating  frequency 

The  carrier  (or  operating)  frequency  of  the  hostile 
signal  is  the  signal  characteristic  of  interest  in  this 
study.  The  intent  of  this  study  is  to  examine  two  methods 
( Blackman-Tukey  and  Burg)  for  determining  the  operating 
frequency  of  one  or  more  emitters,  via  the  power  spectral 
dens ity  ( PSD ) . 
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Background 

The  basic  scenario  is  an  electronic  warfare  (EW) 
receiver  detecting  signals  from  hostile  emitters.  Figure 
1.1  shows  a  generic  block  diagram  of  a  receiver.  No 
distinction  is  made  between  the  receiver  and  the  digital 
processor.  However,  note  that  all  blocks  to  the  right  of 
the  dashed  line  represent  digital  hardware. 


Figure  1.1.  Generic  Receiver  Block  Diagram 

The  receiver  input  signal  r(t)  operates  within  the 
radio  frequency  (RF)  band,  300  MHz  to  30  GHz.  Hence,  the 
signal  is  termed  the  RF  signal.  For  the  purpose  of  this 
study,  the  RF  signal  is  assumed  to  operate  at  the  high  end 
of  the  RF  band  (greater  than  say  IGHz).  The  RF  signal  is 
down-converted  (or  translated)  to  an  appropriate  baseband 
frequency  and  then  sampled,  as  indicated  in  Figure  1.1.  T 
purpose  of  sampling  the  baseband  representation  of  r(t)  is 
to  allow  for  digital  processing.  Generally,  digital 


processing  is  more  efficient  and  less  costly  than  analog 
(28:119).  However,  quantization  error  (or  noise)  is  an 


inherent  problem  of  digital  processing.  This  error  results 
from  the  roundoff  effect  that  occurs  when  converting  a 
continuous  signal  into  a  set  of  finite  discrete  levels. 

Tong  (31:512-521)  provides  a  detailed  discussion  of 
quantization  error  and  methods  for  reducing  this  error. 

The  sampling  theorem  states  that  a  baseband  signal  is 
completely  reconstructed  (no  information  loss)  if  the 
sampling  rate  is  at  least  twice  the  highest  frequency 
component  of  the  baseband  signal  (9:230).  Ideally,  sampling 
is  desirable  at  the  RF  bandpass  signal  -  allowing  for  an  all 
digital  EW  receiver.  Unfortunately,  digital  technology  is 
not  capable  of  providing  reasonable  results  when  the  RF 
signal  is  above  1  GHz.  The  bandpass  sampling  theorem 
requires  that  the  sampling  rate  lie  between  2  BHz  and  4  BHz, 
inclusively.  The  variable  B  is  the  bandwidth  of  the  RF 
signal  which  can  be  on  the  order  of  lOOO's  of  MHz.  Current 
digital  circuitry  is  unable  to  sample  at  such  fast  rates 
(27:55).  Once  the  baseband  signal  is  sampled  and  converted 
into  a  digital  word,  the  digital  computer  is  used  to 
determine  the  PSD  of  the  signal. 

The  above  discussion  did  not  consider  the  statistical 
nature  of  the  bandpass  signal  r(t).  The  signal  r(t) 
consists  of  the  transmitted  hostile  signal  s(t)  plus  noise 
n(t).  The  noise  is  induced  onto  the  hostile  signal  by  the 
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channel  (medium  of  propagation).  Thus,  r(t)  is  actually  a 
random  signal  (or  random  process). 

A  random  process  is  completely  described  statistically 

if  all  its  N-th  order  joint  density  functions  are  known 

(i.e.,  N  approaches  infinity),  a  practical  impossibility 
(9:482).  For  many  practical  applications,  only  the  first 
and  second  order  statistics  are  used  to  describe  a  random 
process  (24:208).  The  autocorrelation  function,  a  second 
order  statistic,  for  a  continuous  time  random  process  is 
defined  as 

R^(t,t+T)  =  E{x(t)x(t+T) }  (1.1 

where  £{•}  denotes  the  expected  value  (or  statistical 
average)  of  the  product  x(t)x(t+T).  If  a  random  process 
x(t)  is  wide  sense  stationary  (WSS),  then 

R^(t,t  +  T)  =  Rj^(t)  (1.2 

This  implies  that  the  autocorrelation  function  is 
independent  of  time  shifts.  The  variable  t  is  called  the 
lag  variable,  the  amount  of  time  by  which  t  lags  behind  t+r 
If  x(t)  is  also  ergodic,  then 

E{ X ( t ) X ( t +T ) }  =  <x(t)x(t+T)>  (1.3 

where  <o>  denote  the  time  average  (13:250).  Thus, 
ergodicity  requires  that  the  result  of  the  statistical  (or 
ensemble)  average  and  the  time  average  be  identical.  The 
receiver  input  signal  r(t)  is  assumed  to  be  both  WSS  and 
ergod i c . 

The  above  considers  the  case  of  a  continuous  time 


random  process.  However,  as  eluded  to  earlier  the 
continuous  process  r(t)  is  down-converted  to  a  baseband 


frequency  and  sampled.  Therefore,  a  discrete  implementation 
of  the  autocorrelation  function  is  desired.  A  sampled 
continuous  time  random  process  is  a  discrete  time  random 
process.  Papoulis  (24:290)  shows  that  if  a  continuous  time 
random  process  x(t)  is  WSS  with  autocorrelation  then 

x(n),  a  sampled  version  of  x(t),  has  an  autocorrelation 
given  by 


R^(k)  =  E{x(n)x(n+k) }  (1.4) 

This  is  called  the  discrete  autocorrelation  function  (or 
autocorrelation  sequence)  and  is  a  sampled  version  of 
R^(t). 

The  ground  work  has  been  developed  and  it  is  now  time 
to  define  the  PSD.  According  to  the  Wiener-Khintchine 
theorem  (14:181),  the  PSD  of  an  ergodic  discrete  time  random 
process  is  given  by 

P^(f )  =  DTFT  (R^(k)  I  (1.5) 

This  says  that  the  PSD  is  the  discrete-time  Fourier 
transform  (DTFT)  of  R^(k)  and,  conversely, 

R^(k)  =  DTFT“‘ip^(  f ) ]  (1.6) 

The  discrete  autocorrelation  function  is  the  inverse 
discrete-time  Fourier  transform  of  the  PSD.  A  more 
definitive  explanation  and  expression  of  Rjj(k)  and  ^re 

presented  in  Chapter  III.  It  is  important  to  understand  the 
motivation  for  stating  the  discrete  autocorrelation 


function.  The  random  signal  processed  by  the  digital 
computer  for  computation  of  its  PSD  is  not  continuous  but 
discrete.  This  is  the  signal  x(n)  that  appears  as  the 
output  of  the  sampler  in  Figure  1.1.  Several  methods  have 
been  developed  for  implementation  on  digital  computers,  that 
determine  the  PSD  of  a  discrete  random  signal  based  on  its 
discrete  autocorrelation  function  (18:1382).  In  order  to  be 
consistent  with  current  literature,  the  discrete  random 
signal  x(n)  is  hence  forth  referred  to  as  a  data  record  (or 
discrete  time  series). 

The  two  methods  considered  in  this  study  are 
Blackman-Tukey  ( BT )  and  Burg.  The  BT  method  is  a  fast 
Fourier  transform  (FFT)  technique  employed  to  determine  the 
PSD  of  a  data  record.  This  method,  commonly  referred  to  as 
a  conventional  approach,  is  computationally  efficient  and 
generally  produces  reasonable  results.  However,  degrading 
frequency  resolution  problems  arise  when  the  data  record 
being  considered  is  short. 

The  Burg  method,  referred  to  as  a  modern  approach, 
attempts  to  overcome  the  frequency  resolution  problems  of 
the  BT  method,  particularly  for  short  data  records.  The 
idea  here  is  to  construct  a  parametric  linear  discrete  model 
that  approximates  the  data  record.  The  PSD  is  computed 
based  on  this  model.  The  motivation  for  using  the  Burg 
method  is  its  promise  for  better  frequency  resolution. 


The  BT  and  Burg  methods  are  discussed  in  detail  in 
Chapter  III. 

Problem  Statement 

Several  approaches  are  available  to  determine  the  PSD 
of  a  random  process.  There  are  pros  and  cons  associated 
with  the  selection  of  any  one  of  the  approaches.  Two 
candidate  approaches  are  examined  in  this  study.  The  two 
approaches  are  applied  to  several  known  data  records.  A 
comparison  of  the  two  approaches  is  made  based  on  their 
statistical  properties  and  computational  efficiency. 

Scope 

This  study  proposes  to  examine  the  merits  of  the  BT  and 
Burg  methods  of  spectral  estimation  for  typical  hostile 
emitter  signals.  The  evaluation  of  the  two  methods  is  based 
solely  on  manual  calculations  and  computer  simulations.  No 
hardware  prototypes  are  attempted.  Also,  a  "smart  routine" 
is  introduced  that  attempts  to  automatically  determine  the 
spectral  peaks  corresponding  to  emitter  frequencies. 

Assumpt ions 

The  emitters  are  assumed  to  transmit  at  a  single 
operating  frequency  per  pulse.  Of  course,  each  emitter 
transmits  at  its  own  respective  frequency.  The  reason  for 
the  restriction  of  a  single  frequency  per  pulse  is  because 
of  software  limitations.  That  is,  the  software  package 
used.  Interactive  Signal  Processing  Executive  (ISPX),  to 


simulate  the  emitter  signals  does  not  allow  the  user  to  vary 
the  frequency  per  pulse.  Thus,  r(t)  is  assumed  to  consist 
of  distinct  frequencies  corresponding  to  each  emitter. 

The  random  process  r(t)  is  assumed  WSS.  Practically, 
many  random  processes  are  not  strictly  WSS.  Hence,  the  PSD 
cannot  be  defined  as  given  by  (1.5).  However,  if  the 
duration  of  the  observation  interval  (the  number  of  data 
points)  is  small,  the  process  may  be  considered  locally  WSS 
(24:125).  A  locally  WSS  process  is  one  for  which  the 
variations  of  the  autocorrelation  function  is  small  over  the 
duration  of  the  observation  interval.  Most  real  world 
random  processes  are  considered  locally  WSS. 

Presentat i on 

Chapter  II  provides  a  brief  discussion  of  the  ISPX 
software  package  used  and  a  discussion  of  previous  research. 
Chapter  III  provides  the  development  of  the  theory  for  both 
the  BT  and  the  Burg  methods  of  spectral  estimation.  Chapter 
IV  discusses  findings  as  they  relate  to  several  known  input 
data  records  for  both  methods.  Chapter  V  provides 
conclusions  and  recommendations  for  further  study. 


II.  ISPX  Software  Package  and  Previous  Research 

I ntroduct i on 

The  purpose  of  this  chapter  is  twofold  -  to  briefly 
discuss  the  ISPX  software  package  and  to  present  previous 
research  on  the  BT  and  Burg  methods  of  spectral  estimation. 
The  ISPX  software  package  currently  resides  on  the  Air  Force 
Wright  Aeronautical  Laboratories  ( AFWAL )  VAX  computer  in  the 
directory  of  Dr.  James  Tsui.  The  package  allows  the  user  to 
generate  a  variety  of  "real  world"  data  records.  A  number 
of  different  spectral  estimation  methods  can  be  applied  to 
each  data  record,  and  the  results  viewed  graphically. 
However,  as  mentioned  in  the  previous  chapter  only  the  BT 
and  Burg  methods  are  considered  in  this  study.  The  ISPX 
package  is  a  tool  which  very  nicely  allows  for  studying  the 
characteristics  of  the  BT  and  Burg  methods.  The  findings 
presented  in  Chapter  IV  are  a  result  of  applying  the 
appropriate  routines  of  ISPX  to  several  known  data  records. 
The  appropriate  routines  are  provided  in  Appendix  A. 
Reference  to  each  routine  is  made  in  the  next  chapter. 

The  data  record  is  traditionally  processed  by  using  FFT 
based  (or  conventional)  methods  of  spectral  estimation. 

There  are  inherent  limitations  associated  with  conventional 
methods.  The  most  prominent  limitations  are  poor  frequency 
resolution  and  leakage.  As  will  be  discussed  in  the 
following  two  chapters,  these  limitations  become  more 
pronounced  for  short  data  records.  Short  data  records  are  a 


conunon  occurrence  for  receiver  systems  like  that  outlined  in 
Chapter  I.  In  an  attempt  to  alleviate  the  inherent 
limitations  of  conventional  methods,  researchers  have 
devised  other  methods  of  spectral  estimation.  The  other 
methods  are  referred  to  in  the  literature  as  modern  spectral 
estimation  methods.  The  modern  methods,  however,  are  not  a 
"fix-all"  solution  to  spectral  estimation  for  all  data 
records.  They  also  have  their  limitations.  The  BT  and  Burg 
methods  are  respectively  considered  conventional  and  modern 
spectral  estimation  methods.  Both  methods  are  discussed  in 
detailed  in  the  next  chapter. 

ISPX  Software  Package 

This  section  is  intended  to  provide  a  brief  overview  of 
the  ISPX  software  package.  The  interested  reader  is 
referred  to  reference  (2:1-20)  for  a  more  complete 
discussion  of  ISPX.  The  ISPX  package  provides  a  number  of 
signal  processing  operations.  However,  only  the  operations 
required  to  analyze  the  previous  two  mentioned  spectral 
estimation  methods  are  discussed.  The  operations  are: 

1.  Generate  data  arrays 

2.  Choose  PSD  method  and  execute 

3 .  Plot  an  array 

The  generate  data  operation  allows  the  user  to  generate 
a  wide  variety  of  "real  world"  data  records.  Each  data 
record  generated  can  consists  of  up  to  five  independent 
sinusoids.  The  variable  parameters  of  each  sinusoid  are 


frequency,  amplitude,  delay,  and  duration.  The  values  of 
these  parameters  are  specified  by  the  user;  however,  default 
values  are  provided.  This  multiple  sinusoid  feature  allows 
the  user  to  simulate  a  data  record  corresponding  to  a 
received  signal  consisting  of  sinusoids  from  different 
emitters.  A  pseudo  random  WGN  process  can  be  embedded  on 
the  data  record  to  emulate  the  effects  of  the  noisy 
channel.  The  noise  amplitude  is  the  parameter  the  user 
changes  to  vary  the  input  signal-to-noise  ratio  (SNR)  of  the 
data  record.  That  is,  for  a  desired  input  SNR,  the  noise 
amplitude  required  may  be  calculated  by  the  following 

— -Ismrn  -j 

where  NA  and  SA  represent  noise  and  signal  amplitudes, 
respectively. 

Once  the  data  record  is  established,  the  PSD  operation 
is  performed.  The  newly  generated  data  record  is  sampled, 
and  the  appropriate  PSD  choice  is  made.  The  user  specifies 
the  sampling  frequency,  such  that  the  Nyquist  baseband 
sampling  theorem  is  obeyed.  A  unity  sampling  frequency 
allows  the  user  to  deal  in  normalized  (or  fractional) 
frequency.  Of  course,  the  frequencies  contained  in  the  data 
record  must  be  chosen  accordingly.  For  simplicity,  the 
examples  presented  in  Chapter  IV  are  all  based  on  a  unity 
sampling  frequency.  ISPX  allows  the  user  to  select  one  of 
several  spectral  estimation  schemes,  both  conventional  and 
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modern.  However,  as  eluded  to  earlier  only  the  conventional 
BT  and  the  modern  Burg  methods  are  considered  in  this  study. 
In  selecting  the  BT  method,  the  user  must  specify  the 
number  of  frequency  samples,  the  largest  lag,  and  the  lag 
window.  Each  of  these  parameters  are  discussed  fully  in  the 
next  chapter.  However,  it  is  worth  noting  the  different  lag 
windows  provided  by  ISPX.  They  are: 

1 .  Cos i ne 

2 .  Hammi ng 

3 .  Kaiser- Bessel 

As  will  be  discussed  in  the  following  chapters,  the  hamming 
lag  window  provides  the  overall  best  performance  for  the 
examples  considered  in  this  study.  In  selecting  the  Burg 
method,  the  user  must  specify  the  number  of  frequency 
samples,  and  the  model  order  number.  Again,  these 
parameters  are  discussed  in  the  next  chapter. 

The  plotting  operation  is  invoked  after  the  user 
specifies  the  appropriate  PSD  method.  The  plotting 
operation  allows  the  user  to  graphically  illustrate  the 
spectral  estimation  schemes.  The  user  is  provided  the 
option  of  plotting  the  results  of  up  to  four  spectral 
estimation  schemes  on  a  single  graph.  This  feature  is 
particularly  useful  for  making  visual  comparisons  between 
different  spectral  estimation  methods,  or  comparisons  of  the 
same  method  with  differing  parameters. 


Previous  Research 

This  section  provides  a  brief  sununary  of  current 
research  in  the  area  of  spectral  estimation.  Specifically, 
research  done  over  the  past  decade  pertaining  tu  the  BT  and 
Burg  methods  of  spectral  estimation  is  presented. 

Cooper  and  Kaveh  (7:313-323)  considers  the  effects  of 
noise  on  the  spectral  estimate  provided  by  the  Burg  method. 
They  note  that  as  the  noise  increases  the  resolution  of  the 
spectral  estimation  becomes  progressively  worse.  In  an 
attempt  to  mitigate  the  effects  of  the  noise,  they  increased 
the  model  order  number  of  the  estimate.  Noting  that  the 
resolution  improved;  however,  spurious  (or  unwanted)  peaks 
began  to  appear.  They  also  compared  the  spectral  estimates 
produced  by  both  the  BT  and  Burg  methods  for  a  given  data 
record.  The  data  record  consisted  of  two  sinusoids 
corresponding  to  doppler  shifts  from  two  radar  targets. 

They  note  that  the  number  of  lags  required  by  the  BT  metfiod 
are  generally  greater  than  the  model  order  number  required 
by  the  Burg  method  to  achieve  a  given  frequency  resolution. 

Haykin  (12:258-262)  applies  the  Burg  method  to  solve 
the  problem  of  estimating  the  angle(s)  of  arrival  of  plane 
wave(s)  impinging  on  a  linear  array  antenna  from  unknown 
d i r ect  i  on ( s  )  .  The  location  of  the  spectral  peaks  indicates 
the  directions  of  the  incident  plane  waves.  For  an  input 
SNR  of  18dB,  he  notes  that  for  a  single  source  (one  target) 
illuminating  a  linear  vertical  array  composed  of  eight 


equally  spaced  horn  antennas:  (1)  the  location  of  the 
spectral  peak  coincides  very  closely  with  the  actual 
direction  of  the  source;  (2)  the  resolution  increases  with 
increasing  model  order  number;  and  (3)  spurious  spectral 
peaks  appear  as  the  model  order  number  is  increased.  In 
fact,  there  may  be  as  many  as  P  peaks  in  the  spectrum.  The 
variable  P  is  called  the  model  order  number  and  will  be 
discussed  in  the  next  chapter. 

Kunt  (19:326-339)  provides  several  examples 
illustrating  the  performance  of  the  BT  spectral  estimator. 

He  notes  that  as  the  data  samples  of  a  given  data  record  are 
decreased,  the  frequency  resolution  of  the  estimator  also 
decreases.  He  investigates  the  resolut i on- leakage 
phenomenon  of  several  lag  windows.  Substantiating  the 
well-known  fact,  that  for  a  sinusoidal  process  the  frequency 
resolution  of  the  estimator  is  largely  determined  by  the 
main  lobe  and  sidelobe  characteristics  of  the  lag  window 
chosen.  These  characteristics  tend  to  become  less  favorable 
as  the  number  of  data  samples  decreases.  For  this  reason, 
the  BT  spectral  estimator  typically  provides  poor  results 
for  short  data  records.  Particularly,  when  the  data  record 
consists  of  two  or  more  closely  spaced  sinusoids. 

Bishop  and  Ulrych  (4:189-192)  investigate  the  Akaike 
criterion  (briefly  discussed  in  the  next  chapter)  for 
determining  the  model  order  number  of  the  Burg  spectral 
estimator.  Ideally,  the  correct  model  number  is  that  number 


which  results  in  a  minimum  value  of  the  Akaike  criterion. 
They  note  that  the  Akaike  criterion  typically  results  in 
model  order  numbers  that  are  too  low.  They  emphasize  that 
model  order  number  selection  is  in  general  a  difficult 
problem,  and  that  the  Akaike  criterion  and  other  proposed 
criteria  serve  only  as  guidelines.  However,  they  state  that 
for  sinusoidal  processes  the  Akaike  criterion  yields  the 
"best"  results . 
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III.  Detailed  Theory 


Introduction 

The  power  spectral  density  (PSD)  and  its  estimation  are 
receiving  much  attention  by  researchers  in  many  disciplines. 
The  PSD  is  considered  an  effective  medium  that  allows  for 
characterizing  a  stationary  random  process  (6:27).  Two 
common  but  distinctly  different  schemes  are  usually 
considered  when  attempting  to  estimate  the  true  PSD.  They 
are: 

1.  The  conventional  method  of  PSD  estimation 

2.  The  modern  method  of  PSD  estimation 

The  most  notable  conventional  method  is  the  method  proposed 
by  Blackman  and  Tukey  ( BT )  -  the  BT  method  (5:187-282). 

This  method  is  based  on  computing  the  DFT  of  a  finite 
windowed  (or  tapered)  autocorrelation  sequence  (ACS).  Thus, 
the  BT  method  produces  a  PSD  estimate  which  is  the  transform 
of  the  window  function  convolved  with  the  simple  PSD 
estimator.  The  window  function  plays  a  very  important  role 
in  determining  the  resolution  and  the  sidelobe  phenomena 
associated  with  the  estimate  of  the  BT  method,  particularly 
for  sinusoidal  processes  embedded  in  additive  WGN .  Often, 
the  selection  of  the  window  function  that  provides  the 
"optimal"  compromise  between  high  resolution  and  low 
sidelobes  is  determined  empirically  (23:84-86).  The  BT 
method  generally  provides  unacceptable  results  for  short 
data  records  -  a  smeared  spectral  estimate  with  very  poor 


resolution  (1:460).  Short  data  records  occur  frequently  in 
practice,  particularly  for  radar  systems  in  which  only  a  few 
data  samples  are  available  for  each  received  pulse 
(18:1381).  A  common  error  is  that  of  appending  a  sequence 
of  zeroes  to  the  data  record  in  hope  for  improved 
resolution.  This  technique,  commonly  referred  to  as  zero 
padding,  merely  reduces  the  spacing  between  adjacent 
spectral  lines  by  interpolation  (21:43). 

The  motivation  for  modern  methods  is  their  promise  for 
better  frequency  resolution,  particularly  for  short  data 
records.  The  most  notable  modern  method  is  the 
autoregressive  ( AR )  method  (11:56).  The  idea  here  is  to 
model  the  data  record  by  an  all-pole  rational  transfer 
function.  This  allows  for  extrapolation  of  data  samples 
beyond  the  observation  interval,  not  implicitly  assuming  the 
samples  are  zero  as  with  conventional  methods  (7:321). 
Theoretically,  the  data  samples  can  be  extended 
indefinitely.  Mote  data  samples  allow  for  computation  of 
more  autocorrelation  lags,  and  therefore  in  general  an 
improvement  in  resolution.  Also,  for  all  practical 
purposes,  the  window  function  is  removed;  therefore,  the 
degrading  effects  of  sidelobes  are  alleviated.  Of  course, 
the  accuracy  of  the  PSD  estimate  provided  by  the  AR  method 
depends  directly  on  the  accuracy  of  the  rational  transfer 
function  used  to  model  the  data  record.  The  approach  for 
estimating  the  PSD,  using  the  AR  method,  is  drastically 
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different  from  that  of  the  BT  method.  The  parameters  of  the 
assumed  model  are  estimated  from  the  available 
autocorrelation  lags,  and  then  the  theoretical  PSD  implied 
by  the  model  is  calculated  using  the  estimated  parameters. 
There  are  several  known  techniques  for  estimating  the 
parameters  of  the  assumed  model  (20:562).  In  this  chapter, 
only  the  Burg  approach  is  considered. 

The  following  sections  provide  a  detailed  discussion  of 
the  BT  and  Burg  methods. 

Conventional  Spectral  Estimation  -  Via  the  BT  Method 

The  PSD  simply  represents  the  power  distribution  of  the 
random  process  along  the  frequency  axis  (26:137).  For  a  WSS 
continuous  random  process  x(t),  the  true  PSD  is  given  by  the 
Fourier  transform 

00 

(  f  )  =  J  (  T  )  exp  (  -  j  2nf  T  )  dr  (3.1) 

-00 

of  the  autocorrelation  function  R^(t).  Therefore,  the  PSD 
of  the  WSS  discrete  random  process  x(n)  is  given  by  the 
d iscrete-time  Fourier  transform  ( DTFT ) 

00 

Pj^(f)  =  ^  Rj^(kT)exp(-j27TfkT)  (3.2) 

k  =  -00 

of  the  ACS  R^(kT),  which  corresponds  to  equal  distant 
shifted  versions  of  (see  Figure  3.1)  (6:291).  The 
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Figure  3.1.  Continuous  and  Discrete  Random  Processes 
and  their  Respective  PSD's  (6:291) 


ACS  {10:54),  assuming  a  real  and  ergodic  x(n),  as  defined 
by  Eg  (1.4)  is  given  as 
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1  im 
M-*qo 


1 

2M  +  1 


M 


^  X  (  n  )  X  ( n  +  1^ ) 
n  =  -M 


(  3 


For  simplicity  in  notation,  the  sampling  interval  T  is 
implied  (i.e.,  k  =  kT  and  n  =  nT).  Properties  of  the  ACS 
are 


R^(k)  = 

R^(0)  >  |R^(k  )  I 


These  properties  indicate  that  the  ACS  is  a  real  even 


sequence  (for  real  x(n))  with  a  maximum  at  the  zero  lag. 


The  computation  of  the  ACS  is  practically  impossible. 


requiring  an  infinite  number  of  data  samples.  An  estimate 


of  the  ACS  is  the  only  alternative.  Hence,  the  PSD  as  given 


by  Eq  (3.2)  can  at  best  be  approximated.  Two  commonly  used 


autocorrelation  estimates  (17:73-74)  are 


N-  |k  1-1 


R^(k)  ^ 


)  x ( n  +  k ) 


(3.5a) 
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R^(k)  = 


1 

N-  Ik  I 


X ( n ) X ( n  +  k  ) 


(3.5b) 


where  N  is  the  total  number  of  data  samples  (i.e.,  x(n)  for 


n  =  0  to  n  =  N-1).  Thus,  a  total  of  2N-1  lags  of  the 


autocorrelation  estimates  are  possible.  The  symbol 


denotes  the  biased  estimator  and  symbol  "  the  unbiased 


estimator.  A  biased  estimator  is  one  in  which  its  expected 


value  (or  mean)  deviates  from  the  actual  value  being 


estimated.  The  expected  value  of  an  unbiased  estimator  is 


the  actual  value  (30:9).  The  expected  values  of  the  biased 


and  unbiased  autocorrelation  estimators  are,  respectively. 


E|i^(k)}  .  [l  - 


(3.6a) 


Rx(k) 


(3.6b) 


.V. 


Observe  that  the  expected  value  of  the  biased 


autocorrelation  estimator  is  equal  to  the  actual  value  ( k j 


weighted  by  the  Bartlett  (or  triangular)  window.  The 


expected  value  of  the  unbiased  autocorrelation  estimator  is 


the  actual  value.  Intuitively,  it  might  be  thought  that  the 


unbiased  autocorrelation  estimator  is  the  best  choice  for 


use  in  Eq  (3.2).  However,  Chen  (6:27-28),  Kunt 


(19:293-295),  and  Marple  (21:146-149)  show  that  as  the  lag 


value  k  approaches  its  limit  (i.e.,  k  approaches  N-1),  the 


variance  (or  fluctuation)  of  the  unbiased  autocorrelation 


estimator  increases  dramatically  as  compared  with  that  of 


the  biased  autocorrelation  estimator.  This  increased 


variance  causes  the  unbiased  autocorrelation  estimator  to 


yield  results  not  always  satisfying  the  properties  of  Eq . 


(3.4).  Figure  3.2  presents  plots  of  both  estimators  for  a 


data  record  generated  by  a  first-order  digital  filter,  with 


a  psuedo  WGN  input  process  (19:  296-297).  The  plots  are  for 


different  values  of  k  and  N.  A  comparison  of  the  plots 


indicates  the  variance  of  the  unbiased  estimator  is  greater 


for  k  ~  N.  However,  the  variance  of  both  estimators  appears 


identical  for  k  «  N.  The  latter  occurs  most  often  in 


typical  radar  type  problems.  Therefore,  some  researchers 


argue  that  the  unbiased  autocorrelation  estimator  is  in  fact 


the  best  choice  based  solely  on  Eq  (3.6). 


A  more  valid  comparison,  suggested  by  Jenkins  and  Watt 


(15:179-181),  is  to  consider  the  mean  square  error  (mse)  of 
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Figure  3.2.  Plots  for  both  Autoc 
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the  two  autocorrelation  estimators.  Jenkins  and  Watt 
concluded  that  the  mse  of  the  biased  autocorrelation 
estimator  ( k )  is  typically  smaller  for  all  values  of  k. 
Therefore,  in  this  study  the  PSD  is  estimated  using  the 
biased  autocorrelation  estimator,  namely 

N-1 

P^(f)  =  ^  R^( k )exp( - j2nfk )  (3.7) 

k  =  -(N-l  ) 

defined  for  -1/2T  <  f  <  1/2T.  This  equation  represents  a 
low-pass  filtered  approximation  of  P^(f)  and  hence  an 
approximation  of  (  f  )  (see  Figure  3.1).  The  use  of  R^  ( k  ) 
in  Eq  (3.7)  is  not  to  suggest  that  the  biased 
autocorrelation  estimator  always  yields  a  superior  estimate 
of  the  PSD.  Jenkins  and  Watt  only  suggested  that  the  PSD 
determined  from  the  biased  autocorrelation  estimator  is 
better  on  average,  particularly  for  large  k. 

The  PSD  estimator  given  in  Eq  (3.7)  is  called  the 
simple  (or  periodogram)  PSD  estimator  (29:238).  This 
estimator,  at  first  glance,  might  appear  to  approach  that  of 
Eq  (3.2)  for  large  N.  If  this  happened,  then  the  simple  PSD 
estimator  is  said  to  be  consistent.  A  consistent  estimator 
is  one  in  which  its  bias  and  variance  both  tend  to  zero  as 
the  number  of  data  samples  increases  (32:71).  Kay 
(17:80-82)  and  Kunt  (19:302-304)  show  that  the  simple  PSD 
estimator  is  not  consistent,  as  a  consequence  of  the 
autocorrelation  estimator.  Although  the  bias  of  the  simple 
PSD  estimator  asymptotically  approaches  zero  as  the  data 


samples  increases,  the  variance 

Var |p^ ( f ) j  S  ( f  )  (3.8 


remains  essentially  constant  -  the  square  of  the  desired 
PSD.  Figure  3.3  shows  the  true  PSD  of  a  random  signal 
consisting  of  a  weak  and  strong  sinusoid  embedded  in  WGN . 


Figure  3.3.  The  Ideal  PSD  of  a  Strong  and  a  Weak  Signal 
Embedded  in  WGN  (19:328) 

Note  that  the  frequency  axis  represents  a  fractional  (or 
normalized)  frequency  -  the  actual  frequency  divided  by  the 
sampling  frequency  f^  =  1/T.  Figure  3.4  depicts  the  result 
of  the  simple  PSD  estimator  for  several  data  record  lengths 
As  expected,  the  varianc<?  appears  to  remain  essentially 
constant  for  different  values  of  N.  Many  researchers  agree 
that  the  simple  PSD  estimator  generally  provides  unreliable 
results,  particularly  if  there  is  no  a  priori  knowledge  of 


the  true  PSD  (1:448). 


A  method  to  reduce  the  variance  of  the  simple  PSD 
estimator  is  functionally  described  in  Figure  3.5.  The 
biased  autocorrelation  estimator  is  computed  from  a  given 
data  record,  windowed,  and  then  transformed. 


Figure  3.5.  A  Functional  Description  of  the  BT  Method 


This  method  is  mathemat ica 1 ly  expressed  as 

M 

^  w(k )R^(k )exp( - j2nfk )  (3.9) 

k  =  -M 

where  M  5  N-1  and  w(k)  is  called  the  lag  window  with  the 
following  properties 

0<w(k)<w(0)=l 

w(k )  =  w( -k )  (3.10) 

w(k)=0  for)k|>M 

The  subscript  BT  is  used  to  designate  this  as  the  Blackman 
and  Tukey  PSD  estimator  -  Blackman  and  Tukey  pioneered  the 
development  of  this  estimator.  The  BT  PSD  estimator  can 
also  be  represented  as 


Pg^(£)  =  DTFTjw(k )R^(k ) 


1/2T 


w(  )p^(er  )d? 


(3.11) 


-1/2T 


where  W(f),  the  spectral  window,  is  the  d i scr e te - t i me 
Fourier  transform  of  the  lag  window.  Thus,  the  BT  PSD 
estimator  is  commonly  considered  a  smoothed  (or  filtered) 
version  of  the  simple  PSD  estimator.  All  is  not  gain, 
however,  when  using  the  BT  PSD  estimator.  The  BT  PSD 
estimator  allows  for  a  reduction  of  variance  at  the  expense 
of  increasing  the  bias,  or  equivalently  a  reduction  of 
sidelobes  at  the  expense  of  deteriorating  the  frequency 
resolution.  Kay  (82-84)  shows  that  the  mean  and  variance  of 
the  BT  PSD  estimator  are,  respectively. 


1/2T 


W(  f-?  )P(^)d^ 


(3.12) 


-1/2T 


Var  PBT(f) 


(3.13) 


k=-M 


If  a  small  bias  is  desired,  then  M  is  chosen  large  such  that 
the  spectral  window  acts  as  a  dirac  delta  function.  On  the 
other  hand,  for  small  variance  M  is  chosen  small  as 
indicated  by  Eg  (3.13).  Clearly,  a  compromise  of  acceptable 


bias  and  variance  has  to  be  determined.  In  most  cases,  this 


compromise  depends  on  the  particular  application  of  the 
estimator.  Figure  3.6  shows  the  BT  PSD  estimator  of  the 
spectrum  given  in  Figure  3.4,  assuming  a  Blackman  lag 


window.  A  comparison  of  Figures  3.6  and  3.6  indicates  that 
although  the  weak  sinusoid  is  not  detected,  the  BT  PSD 
estimator  is  much  smoother  (less  variance).  Therefore,  the 
BT  PSD  estimator  is  considered  more  reliable  (19:334). 

The  selection  of  the  lag  window  is  considered  an  "art" 
in  conventional  PSD  estimation.  The  lag  window  is  not  to  be 
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Figure  3.6.  The  BT  PSD  of  the  Weak  and  Strong  Signal 
Example  (19:334) . 


confused  with  the  data  window.  The  data  window  is 
unavoidable  as  implied  by  Eg  (3.6).  That  is,  the 
autocorrelation  estimator  is  computed  based  on  a  finite 
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number  of  data  samples  (the  data  record).  The 
Blackman-Tukey  approach  as  defined  by  Eq  (3,9)  assumes  that 
the  data  window  is  rectangular.  The  lag  window,  however,  is 
optional  and  is  used  primarily  to  reduce  the  variance 
(11:63).  Unfortunately,  there  is  no  algorithm  that 
precisely  defines  the  approach  for  selecting  the  most 
appropriate  lag  window.  For  this  reason,  many  d iscrete-t ime 
lag  windows  have  been  proposed  for  use  in  PSD  estimation. 

An  empirical  approach  of  window  selection  is  usually  done, 
especially  if  there  is  no  a  priori  knowledge  of  the  random 
process.  It  is  beyond  the  scope  of  this  study  to  consider 
the  detailed  characteristics  of  the  many  proposed  lag 
windows  (see  any  of  the  above  mentioned  references). 

However,  it  suffices  to  say  that  the  spectral  resolution  of 
the  BT  PSD  estimator  depends  largely  on  the  main  lobe  of  the 
spectral  window.  Also,  the  sidelobes  of  the  spectral  window 
greatly  influences  the  variance  of  the  BT  PSD  estimator. 

The  implementation  of  the  BT  PSD  estimator  on  digital 
computers  requires 

M 

Pbt^^I^^  y  w(k)R^(k)exp(-j2rTf  .k  ) 
k  =  -M 

=  DFT  |w(k)R^(k)|  (3.14) 

where  f.=  i/KT  for  0  <  i  <  K-1.  The  discrete  Fourier 
transform  (DFT)  diffe’-s  from  the  DTFT  in  that  both  the  time 
and  frequency  representation  of  the  process  are  discrete. 


29 


Generally,  a  fast  Fourier  transform  (FFT)  algorithm  is  used 
to  compute  Eq  (3.14)  (see  Appendix  A).  The  value  of  K  is 
arbitrary,  but  usually  K  »  M,  so  that  sharp  details  in  the 
PSD  estimate  will  not  be  missed  (21:152).  If  K  »  M,  then 
the  argument  of  Eq  (3.14)  is  zero  padded  from  M+1  to  K. 

Zero  padding  decreases  the  spacing  between  spectral  line 
components  in  the  DFT  and  does  not  really  improve  the 
resolution  between  two  closely  spaced  spectral  components  of 
the  signal  (see  Figure  3.7).  The  spectral  resolution  is 
improved,  particularly  for  short  data  records,  by 
considering  modern  techniques  of  spectral  estimation.  These 
techniques  extend  the  ACS  by  extrapolation  (or  prediction) 
rather  than  appending  zeroes.  The  following  section 
considers  the  modern  technique  of  autoregressive  (AR) 
spectral  estimation. 

Modern  Spectral  Estimation  -  Via  the  AR/Burg  Method 

The  autoregressive  (AR)  spectral  estimation  technique 
is  a  parametric  method  which  attempts  to  model  the  data 
record  with  an  all-pole  rational  transfer  function.  Ulrych 
and  Bishop  (4:185)  argue  that  many  d iscrete-t ime  random 
processes  encountered  in  practice  may  be  described  by 


x(n) 


(3.15) 


'IP'V 


Of 


'jiV  j*  ’•j*''*  •  * 


=  -I 


a  (k)x(n-k)  +  r>(n) 

P 


k  =  l 


where  a  (l),a  (2),. ...a  (P)  are  the  P-th  AR  coefficients, 

p  '  p  '  '  p 


and  r7(n)  is  the  sample  of  a  zero-mean  WGN  random  process 


with  variance  a  .  The  above  difference  equation  is  referred 
P 


to  as  an  AR-process  because  x(n)  is  a  linear  regression  on 
itself  with  r)(n)  representing  the  error.  The  variable  P 
represents  the  order  of  the  AR-process  and  plays  a  very 
significant  role  in  determining  the  spectral  density.  More 
is  said  about  the  significance  of  P  shortly.  Taking  the 
z-transform  on  both  sides  of  Eq  (3.15)  and  combining  terms 
yields  the  system  function 

1 


H(z)  = 


( 3.16  ) 


1  + 


I 


ap(k  )z 


k  =  l 


Thus,  the  AR-process  is  viewed  as  being  generated  by 
applying  a  zero-mean  WGN  process  to  an  all-pole  digital 
filter  (see  Figure  3.8).  The  system  function  is  assumed  to 
be  both  stable  and  causal  (i.e.,  all  poles  lie  inside  the 
unit  circle).  This  condition  is  necessary  to  insure  that 
x(n)  is  WSS  (17:179),  an  assumption  made  throughout  this 
study.  The  evaluation  of  H(z)  along  the  unit  circle,  z  = 
exp(j2rTf),  yields  the  all-pole  transfer  function 

1 


H(f)  = 


(3.17  ) 


1  + 


^  ap( k )exp ( - j 2nf k ) 


k=l 


u  ^  ^  k/a.  L  -'a 
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Figure  3.8.  An  All-Pole  Filter  for  Generating  and  AR 
Processs  (12:250) 


Therefore,  the  PSD  of  the  random  process,  that  implied  by 
the  linear  discrete  model,  is 


(3.18  ) 


^ap(  k)exp(-j2r7fk) 

1^1 


for  -1/2T  <  f  <  1/2T.  This  indicates  that  the  problem  of 

PSD  estimation,  via  the  AR  method,  is  actually  a  problem  of 

parameter  estimation.  The  AR  PSD  is  defined  once  the 

parameters  -la  (1),  a  (2),...,  a  (P),  \  and  P  are 

V  P  P  '  P  Pj 

determined.  A  myriad  of  methods  are  available  for 

estimating  these  parameters.  The  interested  reader  is 

referred  to  Kay  (17:243-261),  Marple  (21:206-230),  Akaike 

(3:716-723),  and  Parzen  (25:723-730).  This  study  only 

examines  the  Burg  method  of  estimating  the  parameters 

/ap(l),  3^(2),...,  ap(P),  ^'pV  Also,  the  procedure  proposed 


by  Akaike  for  estimating  the  model  order  P  is  briefly 


discussed.  However,  before  proceeding  to  discuss  these  two 
methods,  the  basis  for  the  high  frequency  resolution 
attributed  to  the  AR  PSD  method  needs  to  be  made  clear. 


DuBroff  (31:1622-1623)  shows  that  an  equivalent 
representation  of  the  AR  PSD  is  given  by 


k  =  -oo 


(  k  )exp  (  -  j  2rrf  k  ) 


where 


(3.19) 


^^(k) 


R^(k)  for  0  <  k  <  P 

P 

E  a(^)Jt  (k-/)  for  k  >  P 
^=l  ^ 


(3.20) 


Thus,  assuming  that  the  first  P+1  true  autocorrelation  lags 
are  known,  via  Eq  (3.3).  Eq  (3.20)  implies  that 
autocorrelation  lags  may  be  recursively  extended 
indefinitely.  Although,  computing  the  AR  PSD  with  Eq  (3.19) 
is  not  very  practical.  The  usefulness  of  this  equivalent 
form  is  to  readily  illustrate  the  basis  for  the  high 
frequency  resolution.  That  is,  in  contrast  to  the  BT  PSD 
estimator  no  windowing  and  no  zero  padding  occurs.  Hence, 
the  AR  PSD  does  not  possess  the  sidelobe  phenomena  of  the  BT 
PSD  estimator  and  generally  has  much  better  frequency 
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for  developing  the  maximum  entropy  spectral  density. 
Considering  how  best  to  estimate  the  extended 
autocor relat ion  lags.  Burg  argued  that  the  time  series  (or 
data  record)  characterized  by  the  known  and  extrapolatated 
ACS  should  have  maximum  entropy.  This  implies  that  the  time 
series  is  the  most  random.  The  rationale  for  this  choice  is 
that  it  provides  for  the  flattest  and  the  most  minimum  bias 
spectral  estimate  (17:199).  Ulrych  and  Bishop  (4:183-200) 
and  a  host  of  other  researchers  show  that  the  maximum 
entropy  spectral  density  is  identical  to  the  AR  PSD  for 
linear  Gaussian  random  processes.  Hence,  the  terms  AR  PSD 
and  Burg  PSD  are  used  synonymously  in  this  study. 

Burg  also  suggested  that  in  estimating  the  parameters 
fp“>'  . . .  ‘'p}  ,  the  Levinson-Dur bin  recursive 

algorithm  be  used.  A  detailed  discussion  of  this  algorithm 
is  found  in  Kay  (17:181-191).  In  summary,  the 
Levinson-Durbin  recursive  algorithm  provides  estimates  of 
the  AR  parameters  by  the  following  procedure 

N-1 

N  =  0 


For  K  =  1, 2,  .  .  .  ,  P 


and 


ak_i(i)  +  rkak_  (  k-i  )  ,  i  =  l, 

k  ' 


k-1 

(3.22) 


(3.23) 
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The  mathematical  form  of  Eq  (3.23)  is  analogous  to  the 


transmission  of  power  through  a  terminated  two-port  device. 


Hence,  r,  is  known  in  the  literature  as  the  reflection 
k 


coefficient  (12:254) 


^p‘2) . ap(p),  c-p 


1 


The  estimate  AR  parameters  ■{3^(1), 


are  known  once  the  reflection 


coefficient  is  determined.  Burg  proposed  to  estimate  the 
reflection  coefficient  by  minimizing  the  sum  of  the  forward 
and  backward  prediction  error  powers.  The  forward  and 
backward  prediction  errors  are  defined,  respectively,  as 


Fj^(n) 


X  ( n ) 


+ 


^  a^ ( k ) X ( n-k ) 


(3.24a) 


and 


P 

B^(n)  =  x(n-P)  +  ^  ap(k )x(n-P  +  k  )  (3.24b) 

k^l 

These  two  equations  represent  prediction-error  filters 
operating  in  the  forward  and  backward  direction, 
respectively.  The  Burg  method  of  determining  insures 
that  if  operating  in  the  forward  direction,  then  the 
prediction  of  x(n)  given  {x(n-l),  x ( n- 2 ),..., x ( n-P ) }  is  the 
"best"  possible  prediction  (i.e.,the  mean  square  error  or 
forward  prediction  error  power  is  minimized).  Similarly,  if 
operating  in  the  backward  direction,  then  the  prediction  of 
x(n-P)  given  {x(n-P  +  l),  x ( n-P  + 2 ),..., x ( n- 1 )  ,  x(n)}  is  the 
"best"  possible  prediction  (i.e.,  the  backward  prediction 


error  power  is  minimized).  Haykin  (12:243-251)  shows  that 


when  the  forward  and  backward  prediction  error  powers  are 
minimized,  the  resultant  errors,  Eq  (3.24),  take  on  the  form 
of  WGN  processes.  If  this  happens,  then  Eq  (3.24)  is 
equivalent  to  the  AR  process  described  by  Eq  (3.15).  The 
sum  of  the  forward  and  backward  prediction  error  powers  is 
given  by 


P 


k 


+ 


(3.25) 


The  reflection  coefficient  is  obtained  by  minimizing  the 
above  expression.  Before  doing  so,  it  is  convenient  to  make 
the  following  substitution  for  the  forward  and  backward 
error 


F^(n)  =  F^_^(n)  .  r^B^.i(n-l) 


(  3 . 26a  ) 


and 


B^(n)  =  B^.,(n-1)  .  r,F^.i(n) 


(3.26b) 


These  two  relat..ons  provide  a  recursive  method  for 
determining  the  errors  and  are  obtained  by  simply 
.substituting  Eq  (3.22)  into  Eq  (3.24).  After  the 
substitution  of  Eq  (3.26)  into  (3.25),  is  differentiated 
with  respect  to  P^^  and  set  equal  to  zero.  Applying  some 
simple  algebraic  man  i  pu  lat  i  ons ,  it  is  easy  to  show  that  P^^ 
is  given  by 


N-1 

I 


F^_l(n)B^_l(n-l) 


N-1 

.11' 


|2  +  |B^_^(n-l) 


-2 


(3.27) 


where  the  expectations  are  approximated  by  spatial  averages. 


Kay  (17:255-256)  shows  Eq  (3.27)  guaranties  that  <  1. 
Hence,  the  Burg  method  for  determining  the  reflection 
coefficient  provides  poles  which  are  on  or  inside  the  unit 
circle,  resulting  in  a  stable  or  marginally  stable  filter. 

In  an  attempt  to  try  and  put  the  above  method  into  some 
perspective,  the  computations  are  to  proceed  as  follows; 

1.  The  initial  conditions  are: 


O', 


=  R  (0) 

X 


(3.28a) 


and 


FQ(n)  =  BQ(n)  =  x(n),  n  =  0, 1 , 2,  .  .  .  ,N-1  (3.28b) 

2.  Compute  the  reflection  coefficient,  Eq  (3.27),  for 
l<  =  1 

3.  Compute  the  variance,  Eq  (3.23),  for  K  =  1 

4.  Compute  the  AR  coefficients,  Eq  (3.22),  for  )(  =  1 

5.  Compute  the  prediction  error  updates,  Eq  (3.26), 
for  K  =  1 


6.  Increment  K  by  one,  go  bac)<L  to  step  2,  and  repeat. 

7.  Stop  computation  after  the  desired  order  P  is 
reached . 

To  illustrate  these  steps  consider  the  simple  first-order  AR 
process  given  by 

x(n)  =  -aj^(  l)x{n-l)  7?(n)  (3.29) 

Multiplying  both  sides  of  this  equation  by  x(n-lt)  and  talcing 
the  expected  values  for  k  =  0  to  )c  =  l  yields 


R^(l)  R^(0) 


a^(l) 


1 

0 


(3.30) 


This  matrix  notation  is  commonly  referred  to  as  the  set  of 


first-order  Yule-Walker  equations  (4:190).  The  unknown 


parameters  are  obviously  {a^(l),  }.  The  procedure 


outlined  in  the  above  seven  steps  yields  the  following 


relat  ions 


r  =  -2  - 

1  N-1 


X 

N-1  - 

O’""’' 


X  ( n ) X ( n- 1  ) 


(3.31) 


^  +  |x(n-l) 1^ 


^  2  ,  T  I  2  ~  2 

=  (l-KJ  )c>'o 


(3.32) 


a^(l)  = 


(3.33) 


F  ( n )  =  X ( n )  +  r,  x( n-1 ) 

1  L 


(3.34) 


Bj^(n)  =  x(n-l)  +  x(n) 


(3.35) 


Thus,  given  x(n)  for  n  =  0  to  n  =  N-1,  the  estimate 


parameters  are  easily  obtainable.  Note  that 


for  higher  orders  of  P,  steps  two  through  six  are  repeated 


until  the  desired  AR  parameter  estimates  <ap(l),  a^  ,  , 


~  2^4 

a  (P),  r  V  are  obtained.  These  estimates  are  substituted 
P  P  J 


into  Eq  (3.18)  in  order  to  obtain  the  AR  (or  Burg)  PSD 


estimator,  namely 


■  V  ,  '.w  N, 


Unfortunately,  the  recursive  algorithm  outlined  above  does 
not  provide  any  constraint  on  the  order  P.  That  is,  the  AR 
parameter  estimates  may  be  recursively  determined  for  any 
order  P.  However,  if  the  value  of  P  is  too  low,  then  the  AR 
PSD  estimator  results  in  a  smooth  spectral  estimate  with 
poor  frequency  resolution.  On  the  other  hand,  if  the  value 
of  P  is  too  high,  then  the  estimator  results  in  a  spectral 
estimate  with  spurious  peaks.  This  phenomenon  is 
illustrated  in  the  next  chapter. 

A  number  of  researchers,  such  as  Akaike  (3:716-723), 
Parzen  (25:723-730),  and  Kashyap  (32:996-998),  have  proposed 
different  schemes  to  determine  the  "correct"  model  order  P. 
The  word  correct  is  in  quotes,  because  model  order 
determination  is  generally  a  non-trivial  problem  and  in  most 
cases  the  schemes  suggested  serve  only  as  guidelines  (22: 
2-11).  Ulrych  and  Bishop  (4:189-192)  empirically  found 
that  the  procedure  suggested  by  Akaike  provides  the  best 
estimate  of  P  for  AR  linear  Gaussian  random  processes.  This 


procedure  simply  requires  that  P  be  selected  such  that  the 
final  prediction  error  (FPE)  given  by 


r  N  +  p  V  2 

I  N  -  P 


0<P<N/2 


N/2<P<N- 1 


(3.37) 


is  minimized.  The  idea  here  is  to  select  some  maximum  model 
order  L  S  N  -  1.  Then  successively  compute  the  FPE  for 
integer  values  of  P  =  1  to  P  =  L.  The  value  of  P  resulting 
in  the  minimum  FPE  is  the  "correct"  model  order.  This 
approach  is  desirable  if  there  is  no  a  priori  knowledge  of 
the  random  processes.  For  this  study,  however,  model  orders 
are  determined  empirically  since  there  is  such  a  "wealth"  of 
knowledge  about  the  simulated  random  processes  (see  Chapter 


Observe  that  the  denominator  of  Eg  (3.36)  is  simply  the 

squared  magnitude  of  the  DFT  of  the  sequence 

id  (l),...,a  (P)i,  where  a  (0)  =  1.  Thus,  once  the  AR 
\  P  P  J  P 

parameter  estimates  ia  (1),  a  (2),...,  a  (P),  a  are 

I  P  P  P  P  J 

recursively  determined  as  discussed  above,  the  AR  PSD 
estimator,  Eq  (3.36),  is  determined  by  using  a  FFT  algorithm 
( see  Append ix  A) . 

The  statistics  of  the  AR  PSD  estimator  are  generally 
impossible  to  obtain  (17:211).  However,  Makhoul 
(20:568-669)  shows  that  for  large  data  records  (i.e.  N 
approaches  infinity)  the  mean  and  variance  of  the  AR  PSD 
estimator  are,  respectively 


^  Par<'> 


(3.38a) 
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Although  these  results  may  be  of  very  little  practical 
importance.  That  is,  they  may  not  provide  good 
appr ox imat i ons  for  practical  problems.  Nonetheless,  it  is 
instructive  to  note  the  variance  dependence  on  the  model 
order  P.  Observe  that  the  model  order  selection  represents 
for  AR  spectral  estimation  the  classical  trade-off  between 
resolution  and  variance. 

Comparison  of  the  BT  and  Burg  Methods 

A  comparison  of  the  BT  and  Burg  methods  seems  the  logical 
progression  for  the  final  development  of  this  chapter.  Also, 
the  discussion  in  this  section  leads  rather  nicely  to  the 
analysis  presented  in  the  next  chapter. 

Several  areas  are  worth  considering  when  comparing  the 
BT  and  Burg  methods.  The  most  prominent  are:  computational 
complexity,  resolution,  and  variance.  The  Burg  method  is 
generally  more  computationally  burdensome.  This  is  evident 
by  considering  Eq  (3.36).  As  eluded  to  earlier  the 
denominator  in  Eq  (3.36)  can  be  efficiently  evaluated  by  an 
FFT  routine.  The  additional  computation  results  from  having 
to  recursively  determine  the  appropriate  AR  parameter 
estimates  before  involving  the  FFT  routine;  however,  if 
resolution  is  the  only  concern.  Cooper  and  Kaveh  (7:320) 
show  that  the  model  order  requited  by  the  Burg  method  is 


generally  much  less  than  the  number  of  lags  required  by  the 
BT  method.  Under  this  condition,  the  computation  time  of 
the  Burg  method  is  comparable  to  that  of  the  BT  method. 

As  discussed  previously,  the  frequency  resolution  of 
the  BT  method  is  determined  largely  by  the  lag  window  (or 
spectral  window).  A  simple  exercise  in  Fourier  transforming 
readily  shows  that  the  frequency  resolution  in  inversely 
proportional  to  M,  where  M  <  N  -  1  (see  Eq  (3.91).  Marple 
(21:664)  shows  that  the  frequency  resolution  of  the  Burg 
method  for  sinusoidal  processes  is  approximately  given  by 


Af 


AR 


1 .03 _ 

P [ SNR (P  +  1 )  ) ° 


(3.39) 


Thus,  the  frequency  resolution,  via  the  Burg  method, 
decreases  with  decreasing  SNR  but  also  increases  with 
increasing  model  order  number  P.  The  frequency  resolution 
of  the  BT  method  is  totally  independent  of  input  SNR.  More 
is  said  about  frequency  resolution  for  both  methods 
in  the  next  chapter. 

A  comparison  of  Eqs  (3.13)  and  (3.38b)  shows  that  for  a 
given  data  record  length  N  the  model  order  number  P  is 
analogous  to  the  lag  window  length  M.  That  is,  the  variance 
of  both  methods  tends  to  increase  as  their  respective 
parameter  increases. 

The  following  chapter  provides  a  qualitative  comparison 


of  the  two  methods  for  several  data  records. 


IV.  Analysis 


Introduction 

The  purpose  of  this  chapter  is  to  graphically  examine 
the  two  spectral  analysis  techniques  presented  in  the 
previous  chapter.  The  two  techniques  are  applied  to  several 
known  data  records  and  comparisons  made.  Also,  a  "smart 
routine"  is  introduced  that  attempts  to  automatically 
determine  peaks  corresponding  to  actual  frequen^'ies  of 
given  PSD  plots . 

A  major  "limitation"  of  the  receiver  system  described 
in  Chapter  I  is  that  of  pulse  duration.  That  is,  the  pulse 
received  by  the  receiver  is  relatively  short  in  duration, 
resulting  in  only  a  few  data  samples  available  for 
processing.  The  data  record  resulting  from  the  received 
pulse  is,  therefore,  considered  a  short  data  record.  As 
mentioned  in  Chapter  III,  computing  the  PSD  of  short  data 
records  may  provide  less  than  desirable  results. 

The  Avionics  Lab  suggested  that  the  number  64  is  a  good 
average  representative  number  of  data  samples  available  for 
each  received  pulse.  Hence,  each  data  record  considered  in 
this  chapter  consists  of  only  64  data  samples.  Recall,  that 
the  received  pulse  may  be  composed  of  a  number  of  sinusoids. 
That  is,  two  or  more  hostile  emitters  may  be  operating 
simultaneously.  As  stated  in  Chapter  I,  each  emitter  is 
assumed  to  transmit  at  a  single  operating  frequency.  The 
received  pulse,  and  therefore  the  data  record  are  composed 


of  sinusoids  corresponding  to  the  operating  frequency  of  each 
emitter.  For  example,  if  two  emitters  are  operating 
simultaneously  at  different  frequencies,  then  the  data 
record  consists  of  two  equivalent  baseband  sinusoids 
corresponding  to  the  two  emitter  frequencies. 

The  data  records  evaluated  are  intended  to  represent 
"real  world"  type  hostile  emitter  signals.  A  psuedo  random 
WGN  process  is  added  to  each  data  record  and  represents  the 
noise  corruption  caused  by  the  medium  of  propagation. 

Again,  the  BT  and  Burg  methods  of  spectral  estimation  are 
applied  to  each  data  record  and  comparisons  made.  Of  the 
available  lag  windows  provided  by  the  ISPX  software  package, 
the  hamming  window  yields  the  best  compromise  between 
frequency  resolution  and  leakage.  This  window  is  used  in 
determining  the  BT  PSD  estimator  for  all  cases  considered 
with  M  =  N/2  (or  M  =  32),  see  Eq  (3.9).  As  mentioned 
earlier,  determining  the  model  order  number  P  of  an  AR 
process  is  a  non-trivial  task,  and  therefore  will  be 
determined  empirically.  However,  data  records  that  are 
exclusively  sinusoidal  (i.e.,  no  embedded  noise)  require  a 
model  order  number  P  =  2m  for  m  number  of  sinusoids 
(17:213).  This  criterion  will  serve  as  a  basis  for 
empirically  determining  the  correct  model  order  P. 


Problem  One 


The  first  problem  considered  is  a  data  record 
consisting  of  a  single  sinusoid  embedded  in  WGN .  This  input 
data  record  is  given  by 

x(n)=sin[2rT(0.2)n]  +  g{n)  (4.1) 

where  g(n)  is  a  psuedo  random  WGN  process.  Figures  4.1 
through  4.3  show  this  input  process  for  signal-to-noise 
ratios  (SNR's)  of  10,  15,  and  20  dB,  respectively.  Figure 
4.4  shows  the  BT  PSD  estimator  of  Eq  (4.1)  for  each 
specified  SNR.  Observe  that  the  frequency  resolution 
defined  by  the  main  lobe  centered  about  the  fractional 
frequency  0.2  is  independent  of  SNR.  This  point  is  further 
illustrated  by  superimposing  these  three  results  as 
indicated  in  Figure  4.5.  Figure  4.6  depicts  the  result  of 
overlapping  the  Burg  PSD  estimator  for  each  specified  SNR. 
The  model  order  selected  is  two.  Observe  that  as  the  SNR 
decreases  the  respective  spectral  peaks  are  broadened  and 
displaced  from  the  true  position  (indicated  by  the  vertical 
dashed  line).  This  illustration  clearly  shows  that  the 
frequency  resolution  of  the  Burg  PSD  estimator  is  directly 
proportional  to  the  input  SNR.  Furthermore,  it  is  noted 
that  only  the  spectral  peak  corresponding  to  the  20  dB  SNR 
correctly  resolves  the  true  fractional  frequency  of  0.2. 

The  broadening  and  displacement  of  the  other  two  spectral 
peaks  are  due  to  the  increased  noise  (or  lower  SNR).  A 
quantitative  relationship  that  explains  this  phenomenon  is  a 


Relative  Amplitude  (dB) 


SNR^  =  ^15  dB*  Estimator  of  the  Single  Sinusoid  in  WGN 
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Figure  4.4c.  BT  Estimator  of  the  Single  Sinusoid  in  WGN 
SNR  =  20  dB 
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Relative  Amplitude  (dB) 


Iff 


Relative  Amplitude  (dR) 


Figure  4.6.  Superposition  of  the  Burg  Estimator  for  SNR's 
10,  15,  and  20  dB  with  P  =  2 


formidable  task.  In  fact,  very  little  is  known 


quantitatively  about  noisy  AR-processes  (4:184-187). 

However,  it  is  qualitatively  instructive  to  consider 
specific  cases.  As  will  be  shown  in  the  following  examples, 
the  model  order  number  P  can  sometimes  be  increased  to 
mitigate  the  effects  of  noise  for  a  given  SNR.  A  word  of 
caution,  however,  is  in  order  here.  As  discussed  in  the 
previous  chapter  if  the  model  order  number  P  is  too  large, 
spurious  peaks  in  the  spectral  estimate  will  appear.  This 
is  simply  explained  by  considering  Eq  (3.36),  repeated  here 
for  convenience 


P  ,  2 

1  +  ^  a^ ( k ) exp ( - j 2nf k ) 


(4.2) 


Note  that  if  P  is  too  large,  then  extra  poles  are  produced. 
These  extra  poles,  commonly  referred  to  as  noise  poles,  have 
a  tendency  to  situate  themselves  too  close  to  the  unit 
circle  in  the  z-plane,  resulting  in  unwanted  (or  spurious) 
peaks  . 


Problem  Two 

This  problem  and  the  remaining  two  problems  consider 
data  records  with  a  minimum  SNR  of  15  dB.  The  Avionics  Lab 
indicated  that  15  dB  is  typically  the  lowest  SNR  of  interest 
for  their  applications. 

The  problem  addressed  here  is  a  data  record  consisting 
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of  two  unequal  amplitude  sinusoids  embedded  in  WGN .  This 
data  record  is  given  by 

x(n)  =  0  .  2s  i  n  [  2r7(  0 . 2  )  n  ]  +  s  i  n  [  2n  (  0  .  2  5  )  n  ] 

+  q ( n  )  (4.3 

Figure  4.7  shows  this  input  process  for  a  SNR  of  15  dB 
relative  to  the  weaker  sinusoid.  The  BT  PSD  estimator  is 
presented  in  Figure  4.8.  Notice  that  the  main  lobe  of  the 
weaker  signal  is  not  present.  Obviously,  the  BT  PSD 
estimator  is  not  capable  of  distinguishing  both  sinusoids. 
This  is  a  common  dilemma  for  all  conventional  schemes  of 
spectral  estimation,  particularly  for  data  records 
consisting  of  weak  and  strong  sinusoids.  This  is  explained 
by  considering  Eq  (3.9),  repeated  here  for  convenience 

M 

PgT(f>=  ^  w(k  )R^(k  )exp(  -  j2TTfk  )  (4.4 

k  =  -M 

The  BT  PSD  estimator  involves  a  linear  operation  on  a 
weighted  ACS  derived  from  a  given  data  record.  Thus,  the 
spectrum  of  the  sum  of  two  sinusoids  (uncorrelated)  is 
simply  the  sum  of  their  respective  spectra.  The  amplitude 
of  the  spectra  is  directly  proportional  to  the  power  in  the 
sinusoids.  Therefore,  the  amplitude  of  the  spectra 
corresponding  to  the  strong  sinusoid  is  greater  than  that 
corresponding  to  the  weaker  sinusoid.  In  fact,  so  much  so 
that  the  main  lobe  amplitude  of  the  weaker  sinusoid  is 
buried  (or  masked)  by  the  sidelobe  amplitudes  of  the  strong 
sinusoid.  It  is  important,  however,  to  point  out  that  the 
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inability  of  the  estimator  to  detect  both  sinusoids  is  not 
due  to  poor  frequency  resolution.  That  is,  the  two 
sinusoids  are  spaced  sufficiently  in  frequency  to  be 
resolved  by  this  estimator.  Again,  to  reiterate  the  problem 
is  that  the  main  lobe  of  the  weaker  signal  is  canceled  by 
the  sidelobes  of  the  strong  signal.  If  the  problem 
considered  two  equal  amplitude  sinusoids,  then  both 
Sinusoids  would  have  been  resolved.  Recall,  that  the 
resolution  of  the  BT  PSD  estimator  is  govern  largely  by  the 
lag  window  chosen.  Of  course,  a  lag  window  with  a  faster 
roll-off  rate  could  have  been  chosen  possibly  to  retrieve 
the  weaker  sinusoid.  But  as  stated  earlier,  of  the 
available  lag  windows  provided  by  ISPX,  the  hamming  window 
provides  the  best  compromised  between  frequency  resolution 
and  leakage.  Problem  four  considers  an  example  in  which 
sinusoids  are  so  closely  spaced  in  frequency,  that  it  is 
impossible  to  resolve  them  using  the  BT  method  regardless 
of  po we  r . 

The  Burg  PSD  estimator  is  shown  in  Figure  4.9.  As  stated 
earlier,  the  model  order  P  can  be  increased  in  order  to 
mitigate  the  effects  of  noise  for  a  given  SNR.  Observe  that 
for  a  model  order  number  of  four,  the  appropriate  order  for 
no  noise,  the  results  of  the  Burg  PSD  estimator  are  no  better 
than  the  results  of  the  BT  PSD  estimator.  Only  the  strong 
signal  is  detected.  For  model  order  numbers  of  six  and  ten, 
the  Burg  PSD  estimator  is  capable  of  detecting  both  the  strong 
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and  weaker  sinusoids.  Note  the  resolution  of  the 
estimator  is  slightly  improved  for  the  model  order  number 
ten.  However,  a  spurious  peak  is  just  starting  to 
appear.  As  explained  in  the  previous  section,  spurious 
peaks  are  the  result  of  noise  poles.  This  phenomenon  can 
also  be  explained  from  a  statistical  point  of  view  by 
considering  Eg  (3.38b),  repeated  here  for  convenience, 

vat{  a  ^  P^^f)  (4.5) 

It  is  observed  that  as  model  order  P  increases  for  a  given 
data  record  length  N,  the  variance  increases.  Thus,  P  can 
be  considered  analogous  to  the  length  of  the  lag  window  M 
for  the  BT  PSD  estimator,  see  Eqs  (3,12)  and  (3.13). 

Problem  Three 

A  data  record  consisting  of  one  or  more  delayed 
sinusoids  is  another  interesting  and  often  common 
occurrence.  For  simplicity,  a  modified  version  of  the 
previous  example  is  considered.  That  is,  the  strong  signal 
is  delayed  by  sixteen  time  units.  Thus,  the  data  record  is 
given  by 

x(n)  =  0 . 2s  in  (  2rr(  0 . 2  )  n  ]  +  s  i  n  (  2tt  (  0 . 2  5  )  (  n- 16  )  1 

+  g ( n )  (4.6) 

This  process  is  shown  in  Figure  4.10.  Again,  the  BT  PSD 
estimator  is  unable  to  distinguish  both  sinusoids  as 
illustrated  in  Figure  4.11.  A  comparison  of  Figures  4.8  and 
4.11  indicates  that  the  shift  appears  to  have  no  apparent 
effect  on  the  frequency  resolution.  This  should  not  be  too 


surprising  considering  that  the  BT  PSD  estimator  is  a  FFT 
based  operator. 

The  results  of  the  Burg  PSD  estimator  are  shown  in 
Figure  4.12.  Comparing  these  results  with  those  of  Figure 
4.9  indicate  that  the  time  shift  causes  the  required  model 
order  number  to  increase.  A  model  order  number  P  =  20  is 
required  to  resolve  both  sinusoids.  In  the  previous 
example,  the  two  sinusoidal  peaks  first  appeared  at  P  =  6. 
The  reason  for  this  increased  model  order  number  is  due  the 
overall  lower  SNR. 

Problem  Four 

The  final  problem  considers  a  data  record  consisting 
of  four  closely  spaced  sinusoids.  This  data  record  is  given 
by 

x(n)  =  0.2sin[2rT(0.2)nl  +  0 . 7s  in  [  2n  (  0 . 22  )  n  ] 

+  s  in  [  2rT(  0 . 24  )  ]  +  0 . 5s  in  [  2n{  0 . 26  )  ]  +  g(n)  (4.7) 

Figure  4.13  illustrates  this  process.  The  BT  PSD  is  given 

in  Figure  4.14.  Observe  that  the  sinusoids  are  too  closely 

spaced  to  be  resolved  by  this  estimator.  The  main  lobe 

shown  is  a  result  of  the  destructive  and  constructive 

interference  of  the  four  sinusoids  and  in  itself  does  not 

represent  any  one  of  the  sinusoids.  Figure  4.15  presents 

the  Burg  PSD  estimator  for  P  =  24.  Note  that  all  four 

sinusoids  are  resolved  with  no  inaccuracies  cause  by  the 

spurious  peaks.  That  is,  the  highest  spurious  peak  is  about 

24dB  below  the  lowest  desirable  peak.  This  gives  rise  to 
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developing  a  threshold  detecting  scheme  for  determining 
peaks  corresponding  to  actual  frequencies  -  the  focus  of  the 
next  section. 

Threshold  Detection  Routine 

The  above  four  problems  are  by  no  means  an  exhaustive 
list  of  all  possible  signal  types.  However,  the  data 
collected  certainly  suggest  that  the  Burg  method  of  spectral 
estimation  is  far  superior  to  the  BT  method  for  short  data 
records.  This,  of  course,  assumes  that  the  input  SNR  is 
much  greater  than  unity.  Under  the  assumption  that  a 
sufficient  input  SNR  (typically  >  15  dB)  is  available,  the 
method  of  choice  is  obviously  the  Burg  method.  Thus,  the 
threshold  detection  routine  developed  is  based  on  applying 
the  Burg  method  to  short  data  records. 

The  foregoing  analysis  assumes  a  priori  knowledge  of 
the  input  data  records.  This  assumption  is  necessary  in 
order  to  investigate  the  merits  of  both  spectral  estimators. 
In  practice,  however,  a  priori  knowledge  is  generally  not 
known  about  the  input  data  records.  In  this  case,  it  is 
often  very  difficult  to  determine  precisely  the  peaks 
corresponding  to  actual  frequencies  and  the  peaks  resulting 
from  inaccuracies  of  the  spectral  estimator.  For  example, 
given  the  PSD  plot  of  Figure  4.12c  without  having  any 
knowledge  of  the  input  data  record,  it  might  easily  be 
interpreted  as  resulting  from  a  data  record  consisting  of 
six  sinusoids.  Recall,  Figure  4.12c  is  actually  lie  result 


c£  a  data  record  consisting  of  two  sinusoids  at  fractional 
frequencies  of  0.2  and  0.25.  The  other  peaks  are  spurious 
peaks  resulting  from  the  inaccuracies  or  statistics  of  the 
Burg  estimator.  The  obvious  dilemma  is  that  of 
differentiating  between  actual  peaks  and  spurious  peaks. 

One  possible  approach  might  be  that  of  threshold  detection. 

The  problem  proposed  by  the  Avionics  Lab  is  to  devise  a 
routine  that  is  capable  of  automatically  detecting  a  maximum 
of  five  actual  peaks  and  their  corresponding  frequencies. 
That  is,  if  a  data  record  consists  of  m  sinusoids  (unknown 
to  the  observer),  then  is  it  possible  to  detect  up  to  five 
actual  peaks  and  their  frequencies. 

In  order  to  come  up  with  an  "adaptive”  threshold  value, 
several  test  cases  are  investigated  (see  Appendix  B).  The 
test  cases  presented  in  Appendix  B  represent  only  a  few  of 
the  cases  actually  considered.  However,  they  serve  to 
summarize  the  implied  results  of  the  many  cases 
investigated.  The  test  cases  include  data  records,  with  64 
data  samples,  consisting  of  one  to  five  independent 
sinusoids.  The  upper  limit  five  is  a  software  constraint 
(see  Chapter  II).  The  amplitude  (or  power)  of  each  sinusoid 
is  varied.  It  is  observed,  that  if  the  amplitudes  of  the 
sinusoids  vary  by  two  orders  of  magnitude  or  more  (e.g.  a 
data  record  consisting  of  three  sinusoids  with  respective 
amplitudes  of  0.1,  2.0  and  30),  then  retrieving  the  lower 
amplitude  sinusoid  is  virtually  impossible  (see  Figure  B.7). 
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However,  if  the  amplitudes  of  the  sinusoids  vary  by  no  more 
than  one  order  of  magnitude,  then  "all"  sinusoids  may 
possibly  be  retrieved  (see  Figures  B.l  through  B.6).  The 
test  cases  suggest  that  a  suitable  adaptive  threshold  is  28 
dB  below  the  highest  peak  in  the  spectrum.  The  word 
adaptive  is  used  rather  loosely  in  that  a  threshold  of  30  dB 
below  the  highest  peak  does  not  adapt  for  all  cases.  This 
method  of  determining  the  threshold  is  purely  empirical, 
resulting  from  determining  an  average  value  below  the 
highest  peak  for  all  test  cases  considered. 

The  ideal  model  order  number  P  required  to  detect  five 
sinusoids  is  ten.  However,  as  previously  pointed  out  the 
corruption  of  noise  causes  this  number  to  increase.  The 
test  cases  suggest  that  a  model  order  number  of  P  =  24, 
using  the  threshold  value  previously  stated,  allows  the 
estimator  to  correctly  resolve  actual  peaks  in  most  cases 
for  data  records  consisting  of  up  to  five  sinusoids.  It  is 
pointed  out,  however,  that  a  minimum  SNR  of  15  dB  is  used  in 
each  case  considered.  Varying  this  parameter  will  change 
the  required  value  of  P  .  Also  the  software  limitation  of  a 
maximum  of  five  sinusoids  in  a  data  record  certainly  does 
not  represent  all  possible  data  record  types.  In  practice 
it  may  be  necessary  to  change  the  value  of  P  if  it  is 
believed  that  the  actual  data  record  consist  of  many 
sinusoids  (i.e.,  m  »  5).  Obviously,  there  are  several 
parameters  that  effect  the  value  of  P;  thus,  it  is  not 
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possible  to  come  up  with  a  value  suitable  for  all  possible 
scenar ios . 

The  "smart"  routine  that  determines  the  actual  peaks 
and  their  corresponding  frequencies,  incorporating  the  above 
criteria,  is  also  presented  in  Appendix  B.  For  purposes  of 
illustration,  consider  the  PSD  plot  of  Figure  4.16.  Figure 
4.16  is  actually  the  Burg  PSD  of  Problem  Three  for  P  =  24  - 
note  the  many  spurious  peaks.  The  results  of  applying  the 
"smart"  routine  are  as  follows: 
f^  =  0.1992 

and 

f^  =  0.2460 

The  reason  for  the  discrepancy  between  the  computed 
frequencies  and  the  actual  frequencies  is  due  to  way  in 
which  the  Burg  PSD  is  computed.  That  is,  the  Burg  PSD  is 
computed  based  on  256  discrete  data  points  over  the 
frequency  range  0  to  0.5.  Thus,  the  computed  frequencies 
are  actually  multiples  of  0.5/256  (or  1/512)  and  represent 
close  approximations. 

The  following  chapter  presents  conclusions  and 
recommendations  of  the  above  analysis. 


V.  Conclusions  and  Recommendations 

Conclus ions 

This  study  has  investigated  the  performance  of  two 
popular  but  distinctly  different  methods  of  spectral 
estimation  from  an  EW  receiver  point  of  view.  In  an  EW 
environment  the  duration  of  the  received  pulse  is  usually 
relatively  short,  resulting  in  only  a  few  data  samples 
available  for  processing.  Thus,  the  method  used  for 
spectral  estimation  has  to  be  capable  of  providing 
reasonable  results  for  short  data  records.  Several  "real 
world"  data  records,  each  consisting  of  64  data  samples, 
were  analyzed  using  both  the  BT  and  Burg  methods  of  spectral 
estimation.  The  Burg  method  was  found  to  yield  far  superior 
results  in  terms  of  frequency  resolution.  However,  its 
performance  was  determined  to  be  a  function  of  the  input 
SNR.  It  was  noted  that  for  low  input  SNR  (i.e.,  SNR  S  10 
dB)  the  results  of  the  Burg  method  degraded  substantially. 
Another  parameter  affecting  the  resolution  capability  was 
the  model  order  number  P.  Increasing  the  value  of  P  has  a 
tendency  of  increasing  the  frequency  resolution,  as  well  as 
introducing  spurious  peaks  into  the  spectrum.  Therefore, 
the  value  P  has  to  be  selected  carefully.  The  test  cases 
that  were  analyzed  suggest  a  P  =  24  for  detecting  actual 
peaks  of  data  records  consisting  of  up  to  five  independent 
sinusoids.  A  minimum  SNR  of  16  dB  was  used  in  each  test 


case.  Evaluation  of  the  test  cases  suggested  that  an 


appropriate  scheme  for  differentiating  between  actual  peaks 
and  spurious  peaks  was  that  of  threshold  detection. 

Recommendations 

This  study  is  not  intended  to  present  an  all 
encompassing  approach  of  spectral  estimation  for  short  data 
records.  The  two  techniques  presented  in  this  study 
represent  only  a  some  fraction  of  the  many  methods  proposed 
and  currently  being  evaluated  by  several  researchers 
(9:1383).  Several  appropriate  recommendations  are  as 
follows : 

1.  Investigate  some  of  the  other  methods  of 
spectral  estimations.  For  example,  the  method  proposed  by 
Pisarenko  (42:347-366)  which  provides  a  very  accurate 
discrete  spectrum  for  data  records  consisting  of 
deterministic  harmonics  in  WGN . 

2.  Increase  the  maximum  allowable  sinusoids  in  a 
data  record  currently  provided  by  ISPX.  This  will  allow  for 
developing  a  more  comprehensive  data  base;  thus,  allowing 
for  more  conclusive  data. 

3.  Apply  the  two  methods  discussed  in  this  study 
to  other  types  of  data.  For  example,  the  Burg  might  have 
great  promise  if  applied  to  two  dimensional  imaging  data. 

4.  Investigate  and  possibly  incorporate  some  of 
the  schemes  that  have  evoled  of  the  past  decade  for 
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Appendix  A:  Subroutines  Used  for  Spectral  Estimation 

Several  algorithms  are  used  in  the  analysis  sections  of 
this  paper.  This  appendix  provides  a  brief  explanation  of 
each  routine  with  a  FORTRAN  listing.  All  FORTRAN  listings 
appear  exactly  as  they  appear  in  the  ISPX  software  package. 
The  listings  were  not  generated  by  the  author.  They  were 
used  merely  as  a  "tool”  for  analyzing  several  data  records. 


S u b r outine  Blackman-Tukey 

subroutine  black tukey( x, n, mode ,wind,m,nexp,pbt) 

c  This  program  computes  the  Blackman-Tukey  spectral 

c  estimator  as  given  by  (3.9).  Either  the  biased  or 

c  unbiased  autocorrelation  estimator  may  be  used  as  well 
c  as  a  lag  window. 

c  The  spectral  estimate  is  evaluated  at  the  frequencies 
c  F=-l/2+ ( I -1 ) /L  for  1=1,2, ...,L.  The  number  of 
c:  frequencies  is  given  by  L  =  2**NEXP. 

c 

c  Input  Parameters: 

c 

c  x 

c  n 

c  mode 

c 

c  wind 

c 

c  m 

c  nexp 

c 
c 
c 

c  Output  Parameters: 

c 

c  pbt  -Real  array  of  dimension  L=2**NEXPxl  of  samples 

c  of  the  Blackman-Tukey  spectral  estimate,  where 

c  pbt ( i )  corresponds  to  the  spectral  estimate  at 

c  frequency  F=-l/2 + ( I -1 ) /L . 

c 
c 
c 


-Complex  array  of  dimension  Nxl  of  data  points. 

-Number  of  data  points. 

-Set  equal  to  zero  for  unbiased  autocorrelation 
estimator;  otherwise,  biased  estimator  used. 

-Real  array  of  dimension  2M+1  of  lag  window 
weights;  wind(l), . . .,wind(m+l), , . .,wind(2M+l) 
correspond  to  wl -M ],..., w[ 0 ),...,  wl  M ]  . 

-Largest  lag  desired. 

-Power  of  two  which  determines  number  of 
frequency  samples  desired,  L=2**NEXP;  must  be 
chosen  so  that  L  is  >=  2*M+2. 


External  Subroutines: 
PREFFT,FFT 


Notes  : 


The  calling  program  must  dimension  the  arrary 
X , w i nd , pbt . 

The  array  w,r,rcorr,p  must  be  dimensioned  >=  the 
variable  dimension  shown,  or  equal  to  2**NEXP.  Also, 
the  array  RCORR  should  be  dimensioned  >=  2M+1. 

complex  x(l),  w(512),  r(512),  recorr(129) 
dimension  window(l),  p(512),  pbt{l) 
pi  =  4 . *atan ( 1 .  ) 

compute  the  autocorrelation  estimates  from  the  date, 
ml -m+ 1 

call  correlation{n,ml, mode , x,x,rcorr ) 

Window  the  M+1  autocorrelation  estimates  and  insert 
them  into  the  last  M+1  locations  of  rcorr.  Then,  fill 
first  m  points  of  rcorr  array  with  the  complex 
conjugates  of  the  last  m  points  (shift  the 
autocorrelation 

sequence  to  the  right  by  m  samples  so  that  FFT  may  be 
used  )  . 

r(ml )=wind(ml)*rcorr(l) 
do  10  i=l,m 

r (ml+i )=wind(ml+i ) *rcorr( i+1) 
r(ml-i )=wind(ml-i)*conjg(rcorr(i  +  l)  ) 

Zero  pad  the  array  of  windowed  autocorrelation  samples 
to  obtain  an  array  of  dimension  equal  to  L. 

1  =  2  *  *  nexp 

do  20  i=2*m+2, 1=256 
r(i)  =  (0.,0.  ) 

Compute  FFT  of  the  autocorrelation  sequence, 
i nvr  s  =  - 1 
npad  =  1 

call  pr e f ft ( 1, npad, invrs, nexp, w) 
norm=  0 

call  f ft ( 1 , npad , nexp, norm, w, r ) 

Compensate  for  shifting  the  autocorrelation  sequence  to 

the  right  by  m  samples. 

do  30  i=l,l 

f=( i-1 . )/l 

arg=2.*pi*f*m 

p(i)=real(r(i)*cexp( cmplx ( 0 . , ar g ) ) ) 

Transpose  halves  of  FFT  output  so  that  first  PSD  sample 

is  at  a  frequency  of  -1/2. 

do  40  1=1, 1/2 

pbt { i+1/2 ) =p( i ) 

pbt ( i ) =p( i+1/2 ) 

return 


Subroutine  Burg 

subroutine  br ug ( x , n, ip, a , s  i  g 2  ) 

c  This  program  implements  the  Bu*,g  method  for  estimation 
c  of  the  AK  parameters  ( 3 . 2 1 ) - ( 3 . 2 7 ) . 
c 

c  Input  Parameters: 

c 

c  X  -Complex  array  of  dimension  Nxl  of  data  points, 

c  n  -Number  of  data  points, 

c  ip  -AR  model  order  desired, 

c 

c  Output  Parameters: 

c 

c  a  -Complex  array  of  dimension  IPxl  of  AR  filter 

c  parameter  estimates  arranged  as  A(l)  to  A(IP). 

c  sig2  -White  noise  variance  estimate, 

c 

c  Notes : 

c 

c  The  calling  program  must  dimension  the  X,A  arrays, 
c  The  arrays  EFK , EFK 1 , EBK , EBK 1 , AA, RHO  must  be  dimensioned 

c  >=  ( n, n, n, n, ipxip, ip  respectively), 

c 

complex  x(l),a(l),efk(512),ebk(512),efkl(512), 
c  ebk 1 ( 51 2 ) , aa ( 128 , 128 ) , :  sumn,sumd 

dimension  rho(128) 

c  Compute  the  estimate  of  the  autocorrelation  at  lag  zero 

c  (3.21). 

r  ho0  =  0 
do  10  i =1, n 

10  r ho0=r hoO+cabs ( x ( i ) ) **2/n 

c  Initialize  the  forward  and  backward  prediction 

c  er r ors ( 7  .  39  )  . 

do  20  i =2, n 
ef kl ( i  )  =x (  i  ) 

20  ebkl(i-l)=x(i-l) 

c  Begin  recursion, 

do  80  k=l, ip 

c  Compute  the  reflection  coefficient  estimate  (3.27). 

sumn= ( 0 . , 0 . ) 
sumd  =  (  0  . ,  0  .  ) 
do  30  i=k+l,n 

sumn=sumn  +  e  f k 1 (i)*conjg(ebkl(i-l) ) 

30  sumd  =sumd  +  cabs (efkl( i) )**2+cabs(ebkl( i-1) )**2 

aa ( k , k ) =- 2 . *sumn/sumd 

c  Update  the  prediction  error  power  (7.40). 

if(k.eg.l)rho(k)=(l.-cabs(aa(k,k))**2)*rho0 
if(k.gt.l)rho(k)=(l.-cabs(aa(k,k))**2*rho(k-l) 
i  f ( i p . eq . 1 )  go  to  9  0 
if(k.eg.l)  go  to  50 

c  Update  the  prediction  error  filter  coefficients  (3.27). 


40 

c 

50 


do  40  j=l,k-l 

aa(j,k)=aa(j,k-l)+aa(k,k)*conjg(aa(k~j,k-l)) 

Update  the  prediction  error  filter  coefficients  (3.24) 
do  60  i=k+2, n 

efk(i)=efkl(i) +aa (k,k ) *ebkl { i-1 ) 

60  ebk(i-l)=ebkl(i-2)+conjg(aa(k,k)*efkl(i-l) 
do  70  i =k  +  2  ,  n 
efkl(i)=efk(i) 

70  ebkl( i-l)=ebk( i-1) 

80  continue 

c  Find  final  values  of  the  prediction  error  power,  which 
c  is  the  white  noise  variance  estimate,  and  the 
c  prediction  coefficients,  which  are  the  AR  filter 

c  parameter  estimates. 

90  sig2=rho(ip) 
do  100  i  =  l,  ip 
100  a ( i ) =aa ( i , i p ) 
return 
end 


Subroutine  PreFFT 


c 

c 

c 

c 

c 

c 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


subroutine  prefft( n, npad, invr s , nexp, w) 

This  program  sets  up  the  complex  exponential  table 
needed  to  compute  the  fast  Fourier  transform  of  an 
array  of  complex  datasamples  using  a 
dec i mat  ion- in-f requency  algorithm.  Pruning  is 
performed  if  zero  padding  is  requested.  The  output 
table  contained  in  the  array  is  input  to  the  program 
FFT  which  computes  the  fast  Fourier  transform  of  the 
data 

Input  Parameters: 
n  -number  ofdata  samples 

npad  -Set  to  1  for  no  zero  padding  (N-point 

transform),  2  for  double  padding  (2N-point 
transform),  4  for  quadruple  padding  (4N-point 
trans  form)  . 

invrs  -Set  to  -1  for  forward  transform,  1  for  inverse 
trans  form. 

Output  Parameters: 

nexp  -Indicates  power  of  two  exponent  such  that 
n=2**nexp.  Set  to  -1  to  indicate  error 
condition  if  n  is  not  a  power  of  two  in  which 
case  program  terminates  prematurely, 
w  -Complex  array  of  dimension  n*npadxl  containing 

exponential  table. 


c  Notes  : 

c 

c  The  calling  program  must  dimension  the  complex  array  w 

c  greater  than  or  equal  to  n*npad. 

c 

complex  w ( 1 ) , u 
nexp=l 

5  nt=2**nexp 

if  (nt.ge.n)  go  to  10 

nexp=nexp+l 

go  to  5 

10  if(nt.eg.n)  go  to  15 

nexp= -1 
return 

15  nt=n*npad 

ang=8.*atan(l. )/nt 
u  =  cmplx(cos(anq), invrs*sin(ang)  ) 
w(l)  =  (l.,0.  ) 
do  20  i=2,nt 
20  w(i)=w(i-l)*u 

return 
end 


Subroutine  FFT 

subroutine  fft(n,npad,nexp,norm,w,x) 

c  Input  parameters: 

c 

c  n , npad , nexp, w  -See  parameter  list  for  subroutine 

c  "PREFFT" 

c  norm  -Set  to  0  for  forward  transform,  else  the  sum  is 
c  divided  by  n  for  inverse  transform, 

c  X  -Complex  array  of  dimension  nxl  of  data  samples, 

c 

c  Output  parameters: 

c 

c  X  -Complex  array  of  dimension  n*npadxl  of 

c  transform  values, 

c 

c  Notes  : 

c 

c  The  calling  program  must  dimension  arrays  x,w. 

c 

complex  x ( 1 ) , w( 1 ) , t , u 
if(nexp.eq.-l)return 
go  to  ( 30, 20, 10, 10 )npad 
10  n2=n»2 

nt=n2*2 
x(n2  +  l )  =x( 1) 
do  12  K=2,n 
12  x(n2+k)=x{k)*w(k) 


X ( n+ 1 ) =x ( 1 ) 
nx  =  n2  +  1 
X { n+nx ) =x ( nx ) 
j  j  =  3 

do  14  k  =  2 ,  n 

x(n  +  k ) =x(k ) *w( j j  ) 

nx=2n+k 

x(n+nx)=x(nx)*w( jj ) 
j  3  =  j  j  +  2 
mm=  4 

go  to  36 

X ( n  + 1 ) =x ( 1  ) 

do  22  k=2,n 

X ( n+k ) =x ( k ) *w( k ) 

nun=  2 

nt  =  n  *  2 

go  to  35 

nt  =  n 

tnm=l 

ll  =  n 

do  70  k=l,nexp 
nn=ll/2 
j  j=n\m+l 

do  40  i=l, nt, 11 
kk  =  i  +nn 
t  =  X  (  i  )  +  X  (  k  k  ) 
x(kk)=x(i)-x(kk) 

X  (  i  )  =t 

i f ( nn . eq . 1 )  go  to  70 
do  60  j=2,nn 
u  =  w(  j  j  ) 

do  50  i=j,nt,ll 

kk  =  i  +  nn 

t=x ( i ) +x ( kk ) 

x(kk)=(x(i)-x(kk))*u 

x(  i  )=t 

j j j+mm 

ll  =  nn 

mni=nun*  2 

cont i nue 

nv2=nt/2 

nml=nt-l 

j  =  l 

do  90  i  =  l,nn\l 
if(i.ge.j)  go  to  80 
t  =  x(  j  ) 
x( j )=x( i  ) 
x( i ) =t 
k  =  nv2 

if(k.ge.j)  go  to  90 
j  =  j-k 
k  =  k/2 
go  to  85 


f ( norm. eg . 0 )  return 
o  100  i=l,nt 
:(  i  )=x(  i  )/n 
eturn 
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Appendix  B:  Theshold  Detection  Routine 


This  appendix  contains  the  threshold  detection  routine 
called  "Peak  Detector".  The  peaks  are  computed  by 
determining  negative  slope  changes  in  a  given  PSD  plot.  In 
order  to  differentiate  between  actual  peaks  and  spurious 
peaks  an  "adaptive"  threshold  value  was  empirically 
determined  as  discussed  in  Chapter  IV.  Also  several  test 
cases  are  presented 


Routine  Peak  Detector 

Integer  I, K, count 
Character  fname*20 
Complex  z { 512  ) 

Real  Fn,Magn,FnNxt,MagNxt, Slope, SlopeNxt , LMagNxt 
thold, test, ABSFn,PtO,Ptl,MaxPt 
c 

c  The  user  specifies  the  appropriate  filename, 

c 

Print*, 'Enter  a  filename.' 

Read ' ( a ) ' , f name 
c 
c 

Open  ( unit=l, f ile-f name , status= ' old ' ) 
c 

c  Initialization 

c 

count  =  0 
PtO  =  0.0 
Slope  =  0.0 
Fn  =  0.0 
Magn  =  0.0 
c 

c  The  value  of  maximum  peak  is  computed, 

c 

Do  20  K  =  1,257 

Read ( 1, *  )  z ( k ) 

Ptl  =  real ( z ( k ) ) 

MaxPt  =  AMAXl(Ptl,Pt0) 

PtO  =  MaxPt 

LMaxPt  =  10. 0*ALOG10( MaxPt) 


continue 

The  pointer  is  re-initialized 

REWIND  (unit=l) 

The  peaks  are  determined. 

Do  10  I  =  1,257 
FnNxt  =  Fn 
MagNxt  =  Magn 
SlopeNxt  =  Slope 
Readd,*)  Z(I) 

Hagn  =  real (z( I ) ) 

Fn  =  -0.5  +  (I-l)/512.0 
ABSFn  =  ABS(Fn) 

If  (I .EQ.l)  Go  To  10 
Slope  =  (Magn  -  MagnNxt)/(Fn  -  FnNxt) 

If  ( SlopeNxt .GT. 0 . 0 )  then 
test  =  SlopeNxt  +  Slope 

If  (test .LE. SlopeNxt )  then 

LMagNxt  =  10 . 0*ALOG10 ( MagNxt ) 

The  adaptive  threshold  value  is  set. 

thold  =  LMagPt  -  28.0 

If  (LMagNxt .GE. thold )  then 
count  =  count  +  1 

The  actual  peaks  and  their  corresponding 
frequencies  are  printed 

Print*  ABSFn,  LMagNxt,  count 
ENDIF 

ENDIF 

ENDIF 


10 

continue 

1 

c 

CIOSE  (unit=l) 

c 

Stop 

End 

i'y. 


Relative  Amplitude  (dB) 


x(n)  =  0.1sin[2TT(0.1)nl  +  0 . 2s  in  I  2n  (  0 . 1 5  )  n  ] 

+  0 . 3sin[ 2n{ 0 . 2 ) n ]  +  g(n) 


Figure  B.3.  Burg  Estimator  of  Three  Single  Sinusoid  in  WGN 
SNR  =  15  dB  with  P  =  24 
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x(n)  =  0.1sinl2rT{0.1)n]  +  0 . 2s  in  I  2n  (  0 . 15  )  n  ] 

+  0 . 3s  i  n  [  2n  (  0 . 2  )  n  ]  +  0 . 4s  i  n  [  2rT  (  0 . 2  5  )  n  )  +  g(n) 


Figure  B.4.  Burg  Estimator  of  Four  Single  Sinusoid  in  WGN 
SNR  =  15  dB  with  P  =  24 
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0.1sint2n(0.1)nl  +  0 . 2s  i  n  [  2rT  (  0 . 1 5  )  n  ] 
0 . 3sin[ 2n(0. 2 )n 1  +  0 . 4s  in [ 2n( 0 . 25 ) n ) 
0 . 5s  in  I 2n( 0 . 3 ) n ]  +  g(n) 


equency 


Relative  Amplitude  (dB) 


0.1sin[2n(0.1)nl  +  2  .  Os  in  t  2ti  (  0 . 15 )  n  1 
30s  in  I  2tt(  0 . 2  )  n  1  +  g(n) 
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Figure  B.7.  Burg  Estimator  of  Three  Single  Sinusoid  in  WGN 
SNR  =  15  dB  with  P  =  24 
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